home *** CD-ROM | disk | FTP | other *** search
/ The X-Philes (2nd Revision) / The X-Philes Number 1 (1995).iso / xphiles / hp48_1 / least.sq2 < prev    next >
Text File  |  1995-03-23  |  6KB  |  212 lines
  1. Article 5657 of comp.sys.handhelds:
  2. Path: en.ecn.purdue.edu!noose.ecn.purdue.edu!samsung!zaphod.mps.ohio-state.edu!sdd.hp.com!spool.mu.edu!snorkelwacker.mit.edu!bloom-beacon!eru!hagbard!sunic!kuling!henrikj
  3. From: henrikj@kuling.UUCP (Henrik Johansson)
  4. Newsgroups: comp.sys.handhelds
  5. Subject: HP48 Least Squares Approximation
  6. Message-ID: <2014@kuling.UUCP>
  7. Date: 7 Apr 91 09:07:50 GMT
  8. Reply-To: henrikj@kuling.UUCP (Henrik Johansson)
  9. Organization: Upsala University, Sweden
  10. Lines: 198
  11.  
  12.  
  13.       Least Squares Approximations on the HP48
  14.  
  15. In a way the statistic functions on the HP48 are quite good, but
  16. if one wants to do some more complicated curve fitting the built-in
  17. functions are unsatisfactory.  Therefore I sat down and rewrote a
  18. LISP program I wrote a few years ago for more general least squares
  19. approximation problems into HP48 notation.
  20.  
  21. Suppose you have a set of points in the plane and you know that
  22. they will fit well on a line, then you use the built-in functions
  23. to solve the problem.  But if they seem to fit better on a
  24. function like "y = a*x^2 + b*x + c + d*ln(x)" or something?  You don't.
  25.  
  26. This program makes it possible to solve a problem like that above.
  27. The most general form of input 0is
  28.  
  29.   y = a_1 * f_1(x_1, x_2, ..., x_m) +
  30.       a_2 * f_2(x_1, x_2, ..., x_m) +
  31.       ... +
  32.       a_n * f_n(x_1, x_2, ..., x_m)
  33.  
  34. together with a set of data points.  I call the x_i's the independent
  35. variables, and y the dependent variable.
  36.  
  37. Example 1:
  38.  
  39.    x:   1   2   3   4
  40.    y:   3   6  11  18
  41.  
  42. and y is a quadratic function like "a*x^2 + b*x + c" then you enter
  43.  
  44.  Stack level
  45.      4:     [ 1 2 3 4 ]          x values
  46.      3:   [ 3 6 11 18 ]          y values (be sure that they agree by index)
  47.      2:   { 'X^2' X 1 }          f_i's
  48.      1:             'X'          independent variable
  49.  
  50. execute LSQ and get the following result
  51.  
  52.    (symbolic mode on)                       (symbolic mode off)
  53.  Stack level                            Stack level
  54.      2:       .999999999992                 2:       .999999999992
  55.      1: '.99999999967*X^2+                  1: [[ .99999999967 ]
  56.         1.6511E-9*X+                            [ 1.6511E-9 ]
  57.         1.99999999835'                          [ 1.99999999835 ]]
  58.  
  59. The result is that a = 1, b = 0 and c = 2.
  60.  
  61. The number in level 2 tells how well the resulting coefficients (a_i)
  62. make the total function correspond to the given set of data.
  63. A value of 1 is a complete match, and 0 is a complete mismatch.
  64.  
  65. --------
  66.  
  67. It may be the case that you have a `t' that are dependent on more than
  68. one variable.
  69.  
  70. Example 2:
  71.  
  72.    x:     1    -8     5     1     2     3     3
  73.    y:     3     6     2     1     3     0     2
  74.    z:    -5     7     4     1     4     6     1
  75.    t:  -102  -123   108    -2    61    74    10
  76.  
  77. and a reasonable function to try is "t = a*x + b*x*z + c*y + d*y*z + e".
  78. (It is always up to the user to suggest a function!)
  79.  
  80.  Stack level
  81.      4:  [[ 1 3 -5 ]
  82.           [ -8 6 7 ]
  83.           [ 5 2 4 ]
  84.           [ 1 1 1 ]
  85.           [ 2 3 4 ]
  86.           [ 3 0 6 ]
  87.           [ 3 2 1 ]]
  88.      3:  [ -102 -123 108 -2 61 74 10 ]
  89.      2:  { X 'X*Z' Y 'Y*Z' 1 }
  90.      1:  { X Y Z }
  91.  
  92. execute LSQ and get the following result
  93.  
  94.    (symbolic mode on)                       (symbolic mode off)
  95.  Stack level                            Stack level
  96.      2:       .999999999998                 2:       .999999999998
  97.      1: '3.00000000006*X+                   1: [[ 3.0000000006 ]
  98.         3.99999999999*(X*Z)-                    [ 3.9999999999 ]
  99.         6.000000000007*Y+                       [ -6.00000000007 ]
  100.         4.000000000001*(Y*Z)                    [ 4.00000000001 ]
  101.         -7.000000000064'                        [ -7.00000000064 ]]
  102.  
  103. The result is that a = 3, b = 4, c = 6, d = 4 and e = -7.
  104.  
  105. --------
  106.  
  107. As you can see there is two ways to enter the x_i-values.  If there is
  108. only one dependent variable (m = 1) then its values may be entered in
  109. a vector (eg. [ 1 2 3 4 5 ]) but if there are more than one, their values
  110. has to be entered as in the Example 2. The y-values should be entered in
  111. a vector.  If there is only one dependent variable it can be entered
  112. as 'X' instead of { X } in stack level 1.
  113.  
  114. If you want to find a and b in "a*x^b", then you have to do some
  115. pre- and post-processing on the data to make it linear in a and b.
  116.  
  117. It is always the case that `y' is a linear
  118. function of the a_i's (unfortunately).  It had
  119. been nice to have a completely general function...
  120.  
  121. The method I used in the program uses much memory but is somewhat faster
  122. than the other method I could have used (that would have used a lesser
  123. amount of memory).  By this I don't mean that the program is fast...:-)
  124.  
  125. Some parts of the code can surely be rewritten for speed, but I didn't
  126. feel the need to make it better than it is - I'm satisfied with the
  127. result.  Almost.  As can be seen by the examples, the numerical accuracy
  128. could have been better but I don't know how a general method for increasing
  129. the accuracy could be implemented.  If anyone has any ideas - please
  130. let me know.  I haven't consulted any textbooks about this problem so
  131. good advice about what to read would be welcome.
  132.  
  133. If there are any questions, suggestions, comments or bugs - just write
  134. me a mail, and I will answer as soon as I can.
  135. ---------
  136.  
  137. The program takes up 672.5 bytes, and has ckecksum #7543h when stored in a
  138. variable called 'LSQ'.
  139.  
  140. ---- cut here ----
  141. %%HP: T(3)A(R)F(.);
  142. \<<
  143.   \<<
  144.     IF DUP SIZE
  145. SIZE 1 ==
  146.     THEN OBJ\-> 1
  147. SWAP + \->ARRY TRN
  148. CONJ
  149.     END
  150.   \>> \-> fix
  151.   \<<
  152.     IF DUP TYPE 5 \=/
  153.     THEN 1 \->LIST
  154.     END 4 ROLL fix
  155. EVAL DUP SIZE OBJ\->
  156. DROP \-> yv fv vv xv
  157. r c
  158.     \<<
  159.       IF r yv SIZE
  160. 1 GET - c vv SIZE -
  161. OR
  162.       THEN
  163. "Size(s) mismatch"
  164. DOERR
  165.       END 1 fv SIZE
  166.       FOR I xv OBJ\->
  167. DROP fv I GET \-> f
  168.         \<< r 1
  169.           FOR J vv
  170. OBJ\-> DROP 1 c
  171.             FOR K c
  172. 2 + K - ROLL SWAP
  173. STO
  174.             NEXT f
  175. EVAL J 1 - c * 1 +
  176. r + J - ROLLD -1
  177.           STEP
  178.         \>>
  179.       NEXT fv SIZE
  180. r 2 \->LIST \->ARRY DUP
  181. yv fix EVAL DUP 4
  182. ROLLD * DUP 4 ROLLD
  183. SWAP DUP TRN CONJ *
  184. / DUP 4 ROLLD TRN
  185. CONJ SWAP 3 ROLLD
  186. SWAP * OBJ\-> DROP
  187. SWAP DUP 1 CON TRN
  188. SWAP * OBJ\-> DROP SQ
  189. r / SWAP OVER -
  190. SWAP yv DUP DOT -
  191. NEG / SWAP
  192.       IF -3 FC?
  193.       THEN \-> cv
  194.         \<< 0 1 fv
  195. SIZE
  196.           FOR I cv
  197. I 1 2 \->LIST GET fv
  198. I GET * +
  199.           NEXT
  200.         \>>
  201.       END
  202.     \>>
  203.   \>>
  204. \>>
  205. ---- and cut here too ----
  206.  
  207. -------
  208. henrikj@kuling.docs.uu.se               Henrik Johansson
  209. (Nothing fancy about me - at least not in the signature)
  210.  
  211.  
  212.